Stratified spatiotemporal chaos in anisotropic reaction-diffusion systems 
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Numerical simulations of two dimensional pattern formation in an anisotropic bistable reaction-diffusion 
medium reveal a new dynamical state, stratified spatiotemporal chaos, characterized by strong correlations along 
one of the principal axes. Equations that describe the dependence of front motion on the angle illustrate the 
mechanism leading to stratified chaos. 



Pattern formation in nonequilibrium systems has been ex- 
tensively studied in isotropic, two-dimensional media 
Among the most prominent experimental examples are 
Rayleigh-Benard convection and the Belousov-Zhabotinsky 
reaction in the contexts of fluid dynamics and chemical re- 
actions respectively. Recently, there has also been consid- 
erable interest in systems with broken rotational symmetry, 
like convection in liquid crystals |Q] and chemical waves in 
catalytic surface reactions Experimental and theoretical 
studies of such anisotropic systems showed novel phenom- 
ena like ordered arrays of topological defects anisotropic 
phase turbulence [|^], reaction-diffusion waves with sharp cor- 
ners [^], and traveling wave fragments along a preferred ori- 
entation [^. Anisotropy is also often present in pattern forma- 
tion processes in biological media, e.g. in cardiac tissue [^. 

In this Letter we present a new dynamical state that is possi- 
ble only in anisotropic media - stratified spatiotemporal chaos. 
We demonstrate the phenomenon with numerical simulations 
of the bistable FitzHugh-Nagumo equations with anisotropic 
diffusion and characterize it by computing orientation depen- 
dent correlation functions. In addition, an equation for the 
dependence of front velocities on parameters, curvature, and 
the angular orientation is derived. The mechanism leading to 
stratified chaos is described in terms of these analytic results. 
It is tightly linked to the anisotropy of the system and dif- 
fers from the mechanism leading to spiral chaos in isotropic 
bistable media [ll]]. The findings here are relevant to 

catalytic reactions on surfaces where anisotropy is naturally 
provided by crystal symmetry and in biological tissue where 
anisotropy comes from fiber orientation. 

Many qualitative features of pattern formation in chemical 
and biological reaction-diffusion systems are well described 
by FitzHugh-Nagumo (FHN) type models for bistable me- 
dia [ p^ pj[ ]. The specific model we choose to study is 

I ^.-.(„-„»-„) + *-v=« + |[<i*-|l, 



dv ^2 

— = u- aiv - aa + V V , 
at 



(1) 



where u is the activator and v the inhibitor The parameters oi 
and ao are chosen so that Eqs. (|l]) represent a bistable medium 
with two stationary and uniform stable states, an "up" state. 



and a "down" state, V-). Front solutions con- 
nect the two states. The number of front solutions changes, 
when a single front (an "Ising" front) that exists for values of 
rj := ^/eS > rjc loses stability to a pair of counter-propagating 
fronts ("Bloch" fronts) for f] < rjc- The corresponding bifur- 
cation is referred to as the nonequilibrium Ising-Bloch (NIB) 
bifurcation, hereafter the "front bifurcation." The anisotropy 
of the medium is reflected through the parameter d. 

Figure [l| shows the formation and time evolution of strat- 
ified chaos obtained by numerical solution of Eqs. (jlj). The 
initial state is isotropic spiral chaos. As time evolves a clear 
orientation of up-state (grey) domains along the y direction 
develops. The domains consist primarily of elongating stripe 
segments which either merge with other segments or shorten 
by emitting traveling blobs. The stratified chaos state is robust 
and develops from a variety of initial conditions including a 
single spot. 







FIG. 1: The development of stratified spatio-temporal chaos from 
an initial state of isotropic spatio-temporal chaos. The thick contour 
line represents u = and the thin contour line v = 0. The shaded 
regions are up-state domains. Parameters: al — 2.0, ao — —0.1, 
e = 0.039, S = 1.7, d = 1.0, x = [0, 80], y = [0, 80]. 

To gain insight we have followed the dynamics of u and 
V along the x and y axes and display them in the form of 
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space-time plots in Fig. g. A nearly periodic non-propagating 
pattern along the x axis and irregular traveling wave phenom- 
ena along the y axis are observed. A characteristic property 
of spatio-temporal chaotic patterns are correlations that decay 
on a length scale ^ much smaller than the system length L. 
We have computed the normalized spatial two-point correla- 
tion functions, Cy{r) and Cx{r), for the u field in both the x 
and y directions, where Cy{r) —< Au{x, y + t)Au{x, y) > 
I < Au{x, yf >, C^(r) =< Au{x + r, y)Au{x, y) > / < 
Au{x,y)'^ >, Au{x,y) ~ u{x,y)~- < u >, and the brack- 
ets <> denote space and time averaging. Figure || shows the 
results of these computations. Correlations in the y direction 
decay fast to zero, whereas correlations in the x direction os- 
cillate with constant amplitude. This observation may be used 
to define stratified chaos as a state that displays finite correla- 
tion length in one direction (x) and infinite correlation length 
in the other (y). 



anisotropic media as well. 



o 



o 



.X 0.0 




-0.5 



0.0 10.0 20.0 30.0 40.0 50.0 



FIG. 3: Correlation functions Ca;(r') (solid curve) and Cy{r) 
(dashed curve). 




FIG. 2: Space-time plots of cuts parallel to the x axis (top) and 
the y axis (bottom) for the simulation of Fig. |l|. The nearly verti- 
cal columns on the top figure show instances of periodic breathing 
motion of stripe segments (modulated continuous segments) and pe- 
riodic blob formation (spot arrays). 

An important analytical tool for studying front dynamics 
consists of relations between the normal front velocity C,, and 
other front properties like curvature. Relations of that kind 
have successfully been used in the study of pattern formation 
in isotropic media [|lO|, and we wish to exploit this tool for 



Velocity-curvature relations are derived here for A 
\/e/5 ^ 1. The first step in this derivation is to define 
an orthogonal coordinate system (r, s) that moves with the 
front, where r is a coordinate normal to the front and s is 
the arclength. We denote the position vector of the front by 
X(s, t) = (X, Y), and define it to coincide with the u = 
contour line. The unit vectors tangent and normal to the front 
are given by s = cos + sin 6y and f = — sin 6'x + cos 6'y, 
respectively, where ^(s, t) is the angle that s makes with the 
x axis. A point x = (x, y) in the laboratory frame can be ex- 
pressed as X = X(s, t)+rf. This gives the following relations 
between the laboratory coordinates (x, y) and the coordinates 
(s,r) in the moving frame: x = X{s,t) — rsm6(s,t), and 
y = Y{s, t) +r cos 9{s, t) where we defined s = 3X/9s and 
dX/ds = cos 9, dV/ds = sin 6*. In terms of the moving 
frame coordinates the front normal velocity and curvature are 
given by Cn ~ and k — respectively. 

The second step is to express Eqs. (|l]) in the moving frame 
and use singular perturbation theory, exploiting the smallness 
of A. We distinguish between an inner region that includes 
the narrow front structure, and outer regions on both sides of 
the front. In the inner region du/dr ^ ©(A^^) and dv/dr ^ 
C ( 1 ) . In the outer regions both du / dr and dv jdr are of order 
unity. In the inner region v = w/ is taken to be constant. 
Expanding both u and w/ as powers series in A and using these 
expansions in the moving frame equations we obtain, at order 
0(A), a solvabihty condition that leads to the equation 

where = \/l + deos^ 6. In the outer regions to the left 
and to the right of the front region different approximations 
can be made. Here ^ ^ ^ ^ 0(1) and to leading order 
all terms that contain the factor A can be neglected. The re- 
sulting equations can be solved for v in the two outer regions. 
Continuity of v and of ^ at the front position r = yield a 
second relation between C„ and Vf. Eliminating Vf by insert- 
ing this relation into Eq. (^ gives an implicit relation between 
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the normal velocity of the front and its curvature 

^ l + 3/(g)(C„ + ^) ^ 31{0)ao 

(3) 

where = ai + 1/2. A complete account of this derivation 
will be published elsewhere. 

Figures ^ display solutions of Eq. (||) showing the depen- 
dence of the front velocity on front curvature and propagation 
direction for the parameter values of Fig. [l|. In Fig. Qa C„ vs k 
relations are shown for two orthogonal propagation directions. 
In the X direction, 6 — tt/2 (dashed curve), there is only one 
planar front solution with negative velocity (a down-state in- 
vading an up-state), indicating an Ising regime. The negative 
slope of the C„ k relation implies stability to transverse per- 
turbations. In the y direction, 9 = (solid curve), there are 
three planar front solutions indicating a Bloch regime. The 
positive-velocity front (up-state invading down-state) is unsta- 
ble to transverse perturbations whereas the negative-velocity 
front is stable. The middle branch corresponds to an unstable 
front. 
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FIG. 4: Two views of the relation (|3|) corresponding to the parame- 
ters of the simulation in Figure 1 . (a) The velocity-curvature relation 
for fronts in the x direction (dashed curve) and y direction (solid 
curve), (b) The velocity of planar (k = 0) fronts at different angles. 
Bloch fronts exist in narrow sectors around 9 = and 8 = n. The 
wider sectors in between correspond to Ising fronts. 

Figure ^ shows the angular dependence of planar-front ve- 
locities. Counter-propagating (Bloch) fronts exist in a narrow 
sector around the y direction {9 = 0). The other directions 
support only a single (Ising) front. The Ising front speed is 
smallest in the x direction (9 = tt/2) and increases as 9 devi- 
ates from tt/2. 



The dynamics of fronts as displayed in Fig. |l] are affected 
by curvature, propagation direction, and front interactions. 
Eq. (^ captures the effects of the first two factors but does not 
contain information about front interactions. The necessary 
information for our purpose can be summarized as follows. 

The time evolution of a pair of fronts approaching one an- 
other is affected by the speed of non-interacting fronts and by 
the diffusion rate of the activator u (assuming an inhibitor dif- 
fusion constant equal to unity as in Eqs. (|l])). Consider a pair 
of fronts pertaining to up-states invading a down-state, prop- 
agating toward one another in an isotropic medium {d = 0). 
If the distance between the fronts decreases below a critical 
value, Ac ~ 0{^/7JS), the two fronts collapse, leaving a uni- 
form up-state. The accumulation of the inhibitor in the space 
enclosed by the fronts, however, slows their motion. If the 
(non-interacting) front speed is low enough, or if A is small 
enough, there is enough time for the inhibitor to slow the 
fronts down to a complete stop before they reach the criti- 
cal distance Ac. In that case the subsequent evolution depends 
on the front type. A pair of Ising fronts may either form a 
stationary pulse or, closer to the front bifurcation, a breath- 
ing pulse. A pair of Bloch fronts reflect and propagate away 
from one another [pj|]. Thus, high front speed or fast activator 
diffusion {S small and Ac large) lead to collapse, whereas low 
front speed and slow activator diffusion lead to strong front 
repulsion. A similar argument holds for down-state invading 
up-state fronts. 

Returning to anisotropic media we need to know how the 
two factors that affect front interactions, front speed and acti- 
vator diffusion, depend on the direction of propagation. The 
angular dependence of the front speed is already given in 
Fig. The angular dependence of the activator diffusion can 
be deduced by inspecting Eqs. (jl|). It is 1/(5 in the x direction 
and {l + d)/Sin the y direction. Since d > the activator dif- 
fusion constant increases as the propagation direction changes 
from the x direction to the y direction. 

We can discuss now the mechanism of stratified chaos. As 
Fig. Qa shows, the x direction represents a system that sup- 
ports stationary or breathing planar stripe patterns. The y 
direction represents a system where traveling waves prevail. 
The distinct characters of the medium along the two principal 
axes is reflected in the space-time plots of Fig. ^ nearly verti- 
cal columns in the x direction indicate stationary or breathing 
motion, and diagonal stripes in the y direction indicate travel- 
ing wave phenomena. 

The irregular character of the dynamics comes from blob 
formation events as shown in Figure ^. The front speed and 
activator diffusion in the x direction are sufficiently slow for 
Ising fronts to repel one another (rather than collapse). As the 
tip of the stripe segment grows outward and forms a bulge, 
propagation directions deviating from the x direction develop. 
At these directions both the front speed and the activator diffu- 
sion are higher As a result approaching fronts may collapse. 
This is exactly the blob pinching process in frames (c), (d) and 
(e) of Fig. I The process is periodic for extended periods of 
time as indicated by the vertical spot arrays in Fig. g. 

In summary, stratified chaos relies on two main elements: 
(i) stationary or breathing domains vs traveling wave phenom- 



4 





'I 


1 


u 


d y 




^ u 


' u 



FIG. 5: Close up of repeated blob formation. Shaded regions are up- 
state domains. Thick (thin) lines are u — (v = 0) contours. The 
V — line always lags behind the u = line. The tip of a stripe 
segment (a) grows outward (b)-(c). A pinching dynamics begins (d) 
which leads to blob formation (e) traveling along the y direction (f). 
The blob formation leaves a shortened stripe segment (e) whose tip 
grows outward again (f) and the process repeats. The parameters are 
the same as in Fig. hi 



ena in orthogonal directions, and (ii) an angular dependence 
of front interactions that leads to blob formation. Without the 
second element stripe segments would merge to ever longer 
segments until a periodic stripe pattern is formed. These ele- 
ments suggest the parameter regime where stratified chaos is 
expected to be found. The first element implies a regime along 
the front bifurcation such that there is stationary or breath- 
ing motion in the direction with faster activator diffusion and 
traveling waves in the other (slower diffusion) direction. The 
width of this regime increases with the anisotropy d. The sec- 
ond element implies that in the direction of breathing or sta- 
tionary domains the system is close to the onset of breathing 
motion. Deviations from this direction which move the sys- 
tem toward the traveling wave regime then lead to front col- 
lapse. These expectations were verified numerically. Part of 
the work of U. T. was supported by grant D/98/14745 of the 
German Academic Exchange Board (DAAD). M. B. grate- 
fully acknowledges the hospitality of Ben-Gurion University 
and support of the Max-Planck-Society (MPG) by an Otto- 
Hahn-Fellowship. Part of this research is supported by the 
Department of Energy, under contract W-7405-ENG-36. 
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